/******************************************************************************
Safer in School? The Impact of Compulsory Schooling on Maltreatment and Associated Harms
by Adam A. Dzulkipli, Nicole Black, David W. Johnston, and Leonie Segal

Description: This program recreates all of the main tables (i.e. Tables 1 to 4)
presented in the main text of the paper. Tables are exported as .txt files using
the 'esttab' package in Stata. If a table has multiple panels, the panels are
appended together using the 'append' option in 'esttab'. The table will then 
require further manual editing in a typesetting software to get the table to the
exact format as it is presented in the paper. The 'esttab' command can also be 
altered to export the tables as .tex files if so desired.

Note that directory paths have been removed to preserve privacy. Please insert 
the relevant directories in the global macros under "Setup" before running the
.do file. This code was written and implemented in Stata MP 19.5, with a
computational time of approximately <1 minute.

Packages required include:
•	ftools
•	reghdfe
•	estout
•	rdrobust
•	pzms
•	rddisttestk (.ado file by Brigham Frandsen available from
	https://sites.google.com/view/brighamfrandsen/software?authuser=0)
•	rddensity
•	rdhonest

*****************************************************************************/

********************************************************************************
*	Setup
********************************************************************************
clear all
set more off
macro drop all
capture log close

global raw
global tmp 
global clean 
global output 

********************************************************************************
*	Table 1: Covariate balance across the reform cut-off
********************************************************************************

use "${clean}analysis_final.dta", clear

* Construct table in 'wide' format (i.e. estimates presented in columns)
est clear

local h = 721
local rdvars cps_age16_binary sex babywght adv_wgt babygest adv_b_gest adv_cong mother_age adv_m_age adv_m_marstat adv_ses

foreach y in `rdvars'{
	qui eststo: reghdfe `y' post1 h2 c.h2#post1 ///
	if abs(h2) <= `h', vce(robust) absorb(monthofbirth)
	
	qui sum `y' if post1 == 0
	estadd scalar mean = r(mean)
}

* Now we invert the table into 'long' format (i.e. estimates presented in rows)
esttab, keep(post1) noconstant se nostar sca("mean Pre-reform Mean")

matrix C = r(coefs)
matrix S = r(stats)
eststo clear
local rnames : rownames C
local models : coleq C
local models : list uniq models
local i 0
foreach name of local rnames {
    local ++i
    local j 0
    capture matrix drop b
    capture matrix drop se
    foreach model of local models {
        local ++j
        matrix tmp = C[`i', 2*`j'-1]
        if tmp[1,1]<. {
            matrix colnames tmp = `model'
            matrix b = nullmat(b), tmp
            matrix tmp[1,1] = C[`i', 2*`j']
            matrix se = nullmat(se), tmp
        }
    }
    ereturn post b
    quietly estadd matrix se
    eststo `name'
}

local snames : rownames S
local i 0
foreach name of local snames {
    local ++i
    local j 0
    capture matrix drop b
    foreach model of local models {
        local ++j
        matrix tmp = S[`i', `j']
        matrix colnames tmp = `model'
        matrix b = nullmat(b), tmp
    }
    ereturn post b
    eststo `name'
}

* Export table using the 'esttab' package
esttab using "${output}Table1.txt", label replace  ///
	nobase nogaps nomtitle noobs varwidth(40) modelwidth(20) ///
	cells("b(star pattern(1 0 0) fmt(3)) b(star pattern(0 0 1) fmt(3)) b(pattern(0 1 0) fmt(%9.0fc)) se(par pattern(1 0 0) fmt(3))") ///
	star(** 0.05 *** 0.01) ///
	coeflabel(est1 "Any CPS involvement before age 16" est2 "Female" ///
	est3 "Birthweight (in grams)" est4 "Low birthweight (< 2500 g)" ///
	est5 "Gestational age (in weeks)" est6 "Low gestational age (< 37 weeks)" ///
	est7 "Any pediatric conditions" est8 "Mother's age at birth (in years)" ///
	est9 "Mother's age < 21 years at birth" est10 "Mother not married/de facto at birth" ///
	est11 "Low SES (bottom IRSD quintile)") ///
	unstack wide collabels("Coefficient" "Pre-Reform Mean" "Obs." "(Std. Error)") ///
	title("Table 1")
	
********************************************************************************
*	Table 2: RD estimates of the reform's effect on schooling
********************************************************************************

use "${clean}analysis_final.dta", clear

*** Panel A: RD estimates for the entire sample
est clear

local h = 721

foreach y in enrol16 enrol17{
	eststo: reghdfe `y' post1 h2 c.h2#post1 ///
		if rec_census15 == 1 & abs(h2) <= `h', vce(robust) absorb(monthofbirth)
	
	qui sum `y' if rec_census15 == 1 & post1 == 0
	estadd scalar mean = r(mean)
}

* Export table using the 'esttab' package
esttab using "${output}Table2.txt", label replace b(%9.3f) se(%9.3f) noomit ///
	nobase nogaps mtitles("Enrolment at 16" "Enrolment at 17") ///
	keep(post1) ///
	coeflabel(post1 "Reform") ///
	sfmt(3) sca("mean Pre-reform mean") ///	
	star(** 0.05 *** 0.01) ///
	float obslast varwidth(30) modelwidth(20) title("Table 2: Panel A")
	
*** Panel B: RD estimates by past CPS involvement
est clear

local h = 721

foreach y in enrol16 enrol17{	
	eststo: reghdfe `y' nocps_post anycps_post cps_age16_binary h2 c.h2#post ///
		c.h2#cps_age16_binary c.h2#post#cps_age16_binary ///
		if rec_census15 == 1 & abs(h2) <= `h', vce(robust) absorb(monthofbirth#cps_age16_binary)
		
	qui sum `y' if rec_census15 == 1 & post1 == 0 & cps_age16_binary == 0
	estadd scalar mean_nocp = r(mean)
	
	qui sum `y' if rec_census15 == 1 & post1 == 0 & cps_age16_binary == 1
	estadd scalar mean_anycp = r(mean)
}

* Export table using the 'esttab' package (append to Panel A)
esttab using "${output}Table2.txt", label append b(%9.3f) se(%9.3f) noomit ///
	nobase nogaps mtitles("Enrolment at 16" "Enrolment at 17") ///
	keep(nocps_post1 anycps_post1) ///
	coeflabel(nocps_post1 "Reform x No CPS" ///
	anycps_post1 "Reform x Past CPS") ///
	sfmt(3) sca("mean_nocp Pre-reform mean (no CPS)" ///
	"mean_anycp Pre-reform mean (past CPS)") ///	
	star(** 0.05 *** 0.01) ///
	float obslast varwidth(30) modelwidth(20) title("Table 2: Panel B")
	
********************************************************************************
*	Table 3: RD estimates of the reform's effect on new CPS notifications
********************************************************************************

use "${clean}analysis_final.dta", clear

*** Panel A: RD estimates for the entire sample
est clear

local h = 721
local vars cps_after16_binary cps_after16_phy cps_after16_sexual cps_after16_emo cps_after16_neg

foreach y in `vars'{
	qui eststo: reghdfe `y' post1 h2 c.h2#post1 ///
	if abs(h2) <= `h', vce(robust) absorb(monthofbirth)
	
	qui sum `y' if post1 == 0
	estadd scalar mean = r(mean)
}

* Export table using the 'esttab' package
esttab using "${output}Table3.txt", label replace b(%9.3f) se(%9.3f) noomit ///
	nobase nogaps mtitles("All" "Physical" "Sexual" "Emotional" "Neglect") ///
	keep(post1) ///
	coeflabel(post1 "Reform") ///
	sfmt(3) sca("mean Pre-reform mean") ///	
	star(** 0.05 *** 0.01) ///
	float obslast varwidth(30) modelwidth(20) title("Table 3: Panel A")

*** Panel B: RD estimates by past CPS involvement
est clear
scalar drop _all

local h = 721
local vars cps_after16_binary cps_after16_phy cps_after16_sexual cps_after16_emo cps_after16_neg

foreach y in `vars'{
	eststo: reghdfe `y' nocps_post anycps_post cps_age16_binary h2 c.h2#post ///
		c.h2#cps_age16_binary c.h2#post#cps_age16_binary ///
		if abs(h2) <= `h', vce(robust) absorb(monthofbirth#cps_age16_binary)
	
	qui sum `y' if post1 == 0 & cps_age16_binary == 0
	estadd scalar mean_nocp = r(mean)
	
	qui sum `y' if post1 == 0 & cps_age16_binary == 1
	estadd scalar mean_anycp = r(mean)
}

* Export table using the 'esttab' package (append to Panel A)
esttab using "${output}Table3.txt", label append b(%9.3f) se(%9.3f) noomit ///
	nobase nogaps mtitles("All" "Physical" "Sexual" "Emotional" "Neglect") ///
	keep(nocps_post1 anycps_post1) ///
	coeflabel(nocps_post1 "Reform x No CPS" ///
	anycps_post1 "Reform x Past CPS") ///
	sfmt(3) sca("mean_nocp Pre-reform mean (no CPS)" ///
	"mean_anycp Pre-reform mean (past CPS)") ///	
	star(** 0.05 *** 0.01) ///
	float obslast varwidth(30) modelwidth(20) title("Table 3: Panel B")

********************************************************************************
*	Table 4: RD estimates of the reform's effect on emergency department (ED) use
********************************************************************************

use "${clean}analysis_final.dta", clear

*** Panel A: RD estimates for the entire sample
est clear

local h = 721
local vars evered16 evered_17to21 everedinj16 everedment16 evereddigest16 everedoth16
	
foreach y in `vars'{
		qui eststo: reghdfe `y' post1 h2 c.h2#post1 ///
			if abs(h2) <= `h', vce(robust) absorb(monthofbirth)
	
		qui sum `y' if post1 == 0
		estadd scalar mean = r(mean)
}

* Export table using the 'esttab' package
esttab using "${output}Table4.txt", label replace b(%9.3f) se(%9.3f) noomit ///
	nobase nogaps mtitles("At 16" "17 to 21" "Injury" "Ment." "Digestive" "Other") ///
	keep(post1) ///
	coeflabel(post1 "Reform") ///
	sfmt(3) sca("mean Pre-reform mean") ///
	star(** 0.05 *** 0.01) ///
	float obslast varwidth(30) modelwidth(20) title("Table 4: Panel A")
	
*** Panel B: RD estimates by past CPS involvement
est clear

local h = 721
local vars evered16 evered_17to21 everedinj16 everedment16 evereddigest16 everedoth16

foreach y in `vars'{
		eststo: reghdfe `y' nocps_post1 anycps_post1 cps_age16_binary h2 c.h2#post ///
			c.h2#cps_age16_binary c.h2#post#cps_age16_binary ///
			if abs(h2) <= `h' , vce(robust) absorb(monthofbirth#cps_age16_binary)
	
		qui sum `y' if post1 == 0 & cps_age16_binary == 0
		estadd scalar mean_nocp = r(mean)
	
		qui sum `y' if post1 == 0 & cps_age16_binary == 1
		estadd scalar mean_anycp = r(mean)
	}
	
* Export table using the 'esttab' package (append to Panel A)
esttab using "${output}Table4.txt", label append b(%9.3f) se(%9.3f) noomit ///
	nobase nogaps mtitles("At 16" "17 to 21" "Injury" "Ment." "Digestive" "Other") ///
	keep(nocps_post1 anycps_post1) ///
	coeflabel(nocps_post1 "Reform x No CPS" ///
	anycps_post1 "Reform x Past CPS") ///
	sfmt(3) sca("mean_nocp Pre-reform mean (no CPS)" ///
	"mean_anycp Pre-reform mean (past CPS)") ///
	star(** 0.05 *** 0.01) ///
	float obslast varwidth(30) modelwidth(20) title("Table 4: Panel B")	
